#######################################
######Masculinity Threat Analysis######
#######################################
###############10/10/25##################
#######################################

####Analysis Script for Gothreau & Haas, forthcoming in Journal of Experimental Political Science
####A Replication and Extension of Willer et al. (2013) Overdoing Gender: A Test of the Masculine Overcompensation Thesis 

library(haven)
library(ggplot2)
library(dplyr)
library(broom)
library(officer)
library(flextable)
library(forcats)
library(psych)
library(survey)
library(ggplot2)
library(forcats)

#Set Working Directory#
MTData <- read.csv("CleanMTData.csv")

#Create datasets by gender
Data_female <- subset(MTData, gender_binary == "Women")
Data_male <- subset(MTData, gender_binary == "Men")

dvs <- c(
  "IRAQ_WAR_rescaled", "GAY_RIGHT_rescaled", "SUV_DESIRABILITY_rescaled", "SUV_PAY_rescaled",
  "SYSTEM_JUSTIFICATION_rescaled", "DOMINANCE_scale", "TRADITIONALISM_scale",
  "CONSERVATISM_rescaled", "IMMIGRATION_rescaled", "MARIJUANA_rescaled",
  "ELECTRIC_CAR_rescaled", "TRANS_RIGHT_rescaled", "HIRINGPOLICY_WOMEN_rescaled"
)

run_weighted_ttests <- function(data, label) {
  experimental_groups <- setdiff(unique(data$P_EXP), "no threat")
  results <- list()
  
  for (dv in dvs) {
    if (!dv %in% names(data) || !is.numeric(data[[dv]])) {
      warning(paste("Skipping", dv, "- missing or not numeric"))
      next
    }
    
    for (group in experimental_groups) {
      # Subset and clean data
      subset_data <- data %>%
        filter(P_EXP %in% c("no threat", group)) %>%
        select(P_EXP, WEIGHT, all_of(dv)) %>%
        na.omit()
      
      if (nrow(subset_data) < 3) next
      
      # Define survey design
      design <- svydesign(ids = ~1, weights = ~WEIGHT, data = subset_data)
      
      # Run weighted t-test
      formula <- as.formula(paste(dv, "~ P_EXP"))
      ttest <- svyttest(formula, design)
      
      # Extract test statistics
      t_val <- as.numeric(ttest$statistic)
      p_val <- ttest$p.value
      
      # Get weighted means
      group_means <- svyby(as.formula(paste("~", dv)), ~P_EXP, design, svymean, na.rm = TRUE)
      mean_no_threat <- group_means[group_means$P_EXP == "no threat", 2]
      mean_treat     <- group_means[group_means$P_EXP == group, 2]
      
      # Compute difference and estimate standard error
      diff_means <- as.numeric(mean_treat - mean_no_threat)
      se <- ifelse(!is.na(t_val) && t_val != 0, abs(diff_means / t_val), NA)
      
      # Store results
      result <- data.frame(
        gender = label,
        DV = dv,
        experimental_group = group,
        mean_no_threat = round(mean_no_threat, 2),
        mean_treat = round(mean_treat, 2),
        t = round(t_val, 2),
        p = p_val,
        std_error = round(se, 3)
      )
      
      # Significance stars
      result$sig <- ifelse(result$p < .001, "***",
                           ifelse(result$p < .01, "**",
                                  ifelse(result$p < .05, "*", "")))
      
      # Store each comparison result
      results[[paste(label, dv, group, sep = "_vs_")]] <- result
    }
  }
  
  bind_rows(results)
}

ttests_female <- run_weighted_ttests(Data_female, "female")
ttests_male   <- run_weighted_ttests(Data_male, "male")
all_ttests <- bind_rows(ttests_female, ttests_male)

########Treatment Effects#########

# Filter to exclude "threat"
effects_plot_data <- all_ttests %>%
  filter(experimental_group %in% c("constrained", "PCKT"))

# Optional: create human-readable DV labels
effects_plot_data <- effects_plot_data %>%
  mutate(
    DV_label = case_when(
      DV == "IRAQ_WAR_rescaled" ~ "Support for Iraq War",
      DV == "GAY_RIGHT_rescaled" ~ "Negative View of Gay Rights",
      DV == "SUV_DESIRABILITY_rescaled" ~ "SUV Desirability",
      DV == "SUV_PAY_rescaled" ~ "SUV Payment Approval",
      DV == "SYSTEM_JUSTIFICATION_rescaled" ~ "System Justification",
      DV == "DOMINANCE_scale" ~ "Dominance Attitudes",
      DV == "TRADITIONALISM_scale" ~ "Traditionalism",
      DV == "CONSERVATISM_rescaled" ~ "Political Conservatism",
      DV == "IMMIGRATION_rescaled" ~ "Immigration Attitudes",
      DV == "MARIJUANA_rescaled" ~ "Marijuana Legalization",
      DV == "ELECTRIC_CAR_rescaled" ~ "Electric Car",
      DV == "TRANS_RIGHT_rescaled" ~ "Negative View of Trans Rights",
      DV == "HIRINGPOLICY_WOMEN_rescaled" ~ "Gender Affirmative Action",
      TRUE ~ DV
    ),
    effect = mean_treat - mean_no_threat
  )

dodge <- position_dodge(width = 1)

effects_plot_data$gender <- recode(effects_plot_data$gender,
                                   "female" = "Women",
                                   "male" = "Men")

ggplot(effects_plot_data, aes(x = effect, y = fct_rev(DV_label), color = experimental_group)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_errorbarh(aes(xmin = effect - 1.96 * std_error,
                     xmax = effect + 1.96 * std_error),
                 height = 0.2,
                 position = dodge) +
  geom_point(size = 3, position = dodge) +
  facet_wrap(~gender) +
  scale_color_manual(
    name = " ",
    values = c("constrained" = "lightblue", "PCKT" = "darkblue"), 
    labels = c("constrained" = "Constrained Gender Threat",
               "PCKT" = "General Threat")
  ) +
  labs(
    title = " ",
    x = "Treatment Effect (Mean Difference)",
    y = NULL
  ) +
  theme_minimal(base_size = 20) +
  theme(legend.position = "bottom")

ggsave(path = "Visualization", filename = "Mean_byconditions.jpg", width = 15, height = 11, dpi=700)
dev.off()








